Bare Necessities

# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
library(ggh4x)
library(nlme)
library(ggimage)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    axis.line.y = element_line(),
    axis.ticks.y = element_line(), 
    rect = element_rect(fill = "white", color = NA),
    panel.grid.major.y = element_blank(),
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    legend.position = "bottom"))



# labels for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Flag effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * area * season) ) %>% 
    mutate(
      area = factor(group, levels = area_levels_long_all),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, area, season, non_zero))
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

Everything in this markdown uses the size spectra exponent values before I shifted the wmin parameter. The wmin values are fixed at 1g across all groups.

# Data used for Wigley estimates
wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv")) %>% left_join(area_df)

# Load the three regional datasets
wigley_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_bodymass_spectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))

# Load Shelfwide
wigley_bmspectra_shelf <- read_csv(
  here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))



# Join them together
wigley_bmspectra_all <- bind_rows(wigley_bmspectra, wigley_bmspectra_shelf) %>% 
  mutate(metric = "bodymass_spectra",
         community = "Wigley Subset",
         survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  left_join(area_df) %>% 
  mutate(area = factor(area, area_levels_long_all))


# both metrics+communities
spectra_results <-  wigley_bmspectra_all %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels_all),
    metric = fct_rev(metric))

# # Model Datasets for Wigley Species Subset
wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))

# load shelfwide stuff too
wigley_bmspectra_shelf <- read_csv(
  here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))

# Bottom Temperatures
bot_temps <- read_csv(here::here("Data/processed", "trawl_region_seasonal_bottom_temps.csv"))


# Bolt that on
wigley_bmspectra_model_df <- bind_rows(
  wigley_bmspectra_model_df,
  left_join(
    wigley_bmspectra_shelf, 
    bot_temps, 
    by  = join_by("est_year" == "year", survey_area, season))) %>% 
  left_join(area_df) %>% 
  mutate(area = if_else(survey_area == "Northeast Shelf", "Northeast Shelf", area),
         area = factor(area, levels = area_levels_long_all),
  season = fct_rev(season))

# From Andy Beet - Crabs removed
landings_annual <- read_csv(here::here("Data/processed/BEET_GARFO_regional_finfish_landings.csv")) %>% 
  select(survey_area,
         est_year = year,
         total_weight_lb = total_live_lb)

# Get overall landings
shelf_landings_annual <-  mutate(landings_annual, survey_area = "Northeast Shelf") %>% 
    group_by(est_year, survey_area) %>% 
    summarise(total_weight_lb = sum(total_weight_lb, na.rm = T)) %>% 
    mutate(area = survey_area)

Part 1: Summary of Trends

For the first section the focus is on understanding changes in the size structure. This will be done using three metrics:

1. bodymass spectra
   
2. Size quantiles (10th, median, 90th)
   
3. large and small fish indices

This section should address:

Has the size distribution changed over time, where has it changed, and some nuance around what aspects about the distribution have shifted.

For the figures in part 1: non-zero trends over time have been highlighted in the appropriate plot panels with linear trends overlaid.

A. Size Spectra Trends

Size spectra are steepening in 4 seasonal communities. Gulf of Maine (Spring and Fall), Georges Bank (Spring), and Southern New England (Spring).

The size distribution in the Mid-Atlantic Bight (Spring) is becoming more shallow.

# Model linear trends in time
spectra_trend_mod <- gls(
  b ~ scale(est_year) * area * season,
  correlation = corAR1(form = ~ est_year | area/season),
  data = wigley_bmspectra_all)

#check_model(spectra_trend_mod)
# plot(check_collinearity(spectra_trend_mod))


# Get the predictions and flag whether they are significant
bmspectra_trend_predictions <- get_preds_and_trendsignif(spectra_trend_mod) %>% 
  mutate(model_id = "bodymass_spectra-Wigley Subset") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric),
         area = factor(area, levels = area_levels_long_all))


# Make the plot
bmspectra_trend_p <- ggplot() +
   geom_ribbon(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = wigley_bmspectra_all,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = wigley_bmspectra_all,
    aes(est_year, b, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries" ) +
  geom_line(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(area ~ ., scales = "free") +
  labs(
    x = "Year",
    y = "Mass Distribution Exponent",
    linetype = "Trend",
    color = "Season",
    fill = "Season")
 
# Show plot
bmspectra_trend_p

B. Median Weight Trends

Median body weight is declining in 3 seasonal communities, Gulf of Maine (Spring and Fall) and Georges Bank (Spring).

Median weight in the Mid-Atlantic Bight has been increasing in the Spring.

# Load the median weight/length data
wigley_size_df <- read_csv(here::here("Data/processed/wigley_species_size_summary.csv"))

# shelf-wide summary
wigley_sizes_shelf <- read_csv(here::here("Data/processed/shelfwide_wigley_species_size_summary.csv"))


# Join them
size_results <- bind_rows(wigley_size_df, wigley_sizes_shelf) %>% 
  mutate(
    community = "Wigley Subset",
    survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  left_join(area_df)



# Get trends:
wt_mod_wig <- gls(
  med_wt_kg ~ scale(est_year) * area * season,
  correlation = corAR1(form = ~ est_year | area/season),
  data = drop_na(size_results, med_wt_kg))



# Get the predictions and flag whether they are significant
size_trend_predictions <- get_preds_and_trendsignif(wt_mod_wig) %>% 
  mutate(model_id = "med_wt_kg-Wigley Subset") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))




# Median length - all species
size_long <- size_results %>% 
  pivot_longer(cols = c(med_len_cm, med_wt_kg), 
               names_to = "metric", 
               values_to = "val")  %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels_all),
    area = fct_relevel(area, area_levels_long_all),
    season = fct_rev(season))


# Median weight - wigley species

# weight plots
wt_plot <- size_long %>% 
  filter(metric == "med_wt_kg") %>% 
  drop_na(val) %>% 
  ggplot() +
  geom_ribbon(
    data = filter(size_trend_predictions, metric == "med_wt_kg", non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, val, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, val, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  geom_line(
    data = filter(size_trend_predictions, metric == "med_wt_kg", non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(area ~ .,
             scales = "free")  +
  labs(
    x = "Year",
    y = "Median Weight (kg)",
    color = "Season",
    linetype = "Trend",
    fill = "Season")

wt_plot

Notes:

There are a number of years when median weight surges nearly 0.5kg in the Mid-Atlantic Bight. For these situations I think we likely could trace these down to exceptionally high catches of larger species like elasmobranchs.

We are not seeing a see a shelf-wide decline in size so I might need to check what Friedland did or prepare for a confrontation.

C. Bottom Temperature Trends

Bottom temperature data is from the DuPontavice and Saba combined bottom temperature product. This data product has been clipped to the study areas and seasonally averaged to match the survey sampling seasons.

Plotted below are the long-term trends in each area for each survey season. Significant seasonal trends are shown with linear trends and error displayed.

Seasonal bottom temperatures have been increasing in all 8/10 seasonal communities.

Spring temperature in SNE + MAB have not increased over this period.

# Model linear trends in time
btemp_mods <- lm(
  formula = bot_temp ~ scale(est_year) * area * season,
  data = wigley_bmspectra_model_df)



# change the stupid labels around because
btemp_preds <- as.data.frame(
    # Model predictions
    ggpredict(
        btemp_mods, 
        terms = ~ est_year * area * season) ) %>% 
    mutate(
        area = factor(group, levels = area_levels_long_all),
        season = factor(facet, levels = c("Spring", "Fall")))
btemp_signif <- emtrends(
      object = btemp_mods,
      specs =  ~ area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
btemp_trend_predictions <- left_join(btemp_preds, btemp_signif)  %>% 
  mutate(area = group,
         area = factor(area, levels = area_levels_long_all),
         season = factor(season, levels = c("Spring", "Fall")))




# Plot temperatures
temp1 <- wigley_bmspectra_model_df %>% 
  ggplot() +
  geom_ribbon(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, bot_temp, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, bot_temp, color = season, linetype = "5-Year Moving Average"),
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  geom_line(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_color_gmri() + 
  scale_fill_gmri() + 
  #scale_linetype_manual(values = c(1,1)) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  labs(
    y = "Bottom Temperature", 
    x = "Year", 
    color = "Season",
    fill = "Season", 
    linetype = "Trend",
    title = "")



temp1 

Seasonal bottom temperatures have increased in both spring and fall for the Northeast Shelf.

Looking at different areas within the region, Fall temperatures are increasing across the board, but seasonally spring temperatures in SNE and MAB are stable.

D. All Trends Figure

# Show plot
( (wt_plot + theme(strip.text.y = element_blank(),
   plot.margin = margin(0,2,0,2))) | 
    (bmspectra_trend_p + theme(strip.text.y = element_blank(),
   plot.margin = margin(0,2,0,2)))  |
    temp1 + theme(plot.margin = margin(0,2,0,2)) ) +
  plot_layout(guides = "collect")

E. Small/Large Fish Indices

NOTE: Percentile location of the size distribution are determined using DescTools::Quantile() which uses numlen_adj to weight everything by abundance, properly accounting for abundance.

The following figure shows changes in seasonal large and small fish indices within each region. The large and small size thresholds are set for each species independently using their length distributions as in the table above.

The axes on these are: \[LFI = \frac{\sum({N_{tjsk} > size~threshold_k})} {\sum{N_{tjs}}}\]

for each: \(year_t\), \(region_j\), \(season_s\) & \(species_k\).

# Get percentiles based on length for each species
# # Get 90th or 95th percentile as size threshold
# Do them separately
species_size_quants <- wigley_in %>% 
  mutate(area = "Northeast Shelf") %>% 
  unite("species_area", c(area, comname), sep = "X", ) %>% 
  split(.$species_area) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "species_area") %>% 
  separate(col = species_area, into = c("area", "comname"), sep = "X") %>% 
  mutate(area = factor(area, levels = area_levels_long_all)) %>% 
  arrange(area)


# # Show the table
# species_size_quants %>% 
#   mutate_if(is.numeric, round, 2) %>% 
#   head(n = 10) %>% 
#   gt() %>% 
#   tab_header(title = "Wigley Community - Body Mass Quantiles (kg)") %>% tab_spanner(label = "Percentile", columns = -area) %>% 
#   cols_label("area" = "")

# Assign large/small indicators
wigley_lfi <- wigley_in %>% 
  left_join(select(species_size_quants, comname, `10%`, `90%`)) %>%   
  mutate(
    small = ifelse(ind_weight_kg < `10%`, T, F),
    large = ifelse(ind_weight_kg > `90%`, T, F))


# Get total abundance in each area*year*season
total_abundbio <- bind_rows(
    mutate(wigley_lfi, area = "Northeast Shelf"),
    wigley_lfi) %>% 
  group_by(area, est_year, season) %>% 
  summarise(
    total_abund = sum(numlen_adj, na.rm = T),
    .groups = "drop")
  
  
# Get large/small fish  abundances
lf_abundbio <- bind_rows(
    mutate(wigley_lfi, area = "Northeast Shelf"),
    wigley_lfi) %>%
  group_by(area, comname, est_year, season) %>%
  summarise(
    small_abund = sum(numlen_adj * small),
    large_abund = sum(numlen_adj * large),
    .groups = "drop") %>% 
  group_by(area, est_year, season) %>% 
  summarise(
    small_abund = sum(small_abund, na.rm = T),
    large_abund = sum(large_abund, na.rm = T),
    .groups = "drop")

# Join and divide
lfi <- left_join(total_abundbio, lf_abundbio) %>%
  mutate(
    LFI = large_abund / total_abund,
    SFI = small_abund / total_abund)  %>%
  mutate(
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) 

# Plot them as loing-term average

# Plot the SFI
sfi_p <- lfi %>% 
  ggplot() +
  geom_point(
    aes(est_year, SFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, SFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  theme(strip.text.y = element_blank()) +
  facet_grid(area~., scales = "fixed", labeller = label_wrap_gen(width = 8)) +
  labs(
    x = "Year",
    y = "Percent of Total Abundance",
    title = "Wigley Species",
    subtitle = "Small Fish Index",
    fill = "Season",
    color = "Season",
    linetype = "")


# Plot the LFI
lfi_p <- ggplot(lfi) +
  geom_point(
    aes(est_year, LFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, LFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(
    labels = label_percent(), 
    expand = expansion(add = c(0,0.05))) +
  facet_grid(area~., scales = "fixed", labeller = label_wrap_gen(width = 8)) +
  labs(
    x = "Year",
    y = "",
    subtitle = "Large Fish Index",
    fill = "Season",
    color = "Season",
    linetype = "")



# Plot together
(sfi_p | lfi_p) + plot_layout(guides = "collect")

From these figures we can see that there is a general pattern where larger individuals (relative to their conspecific size-distribution) make up a smaller fraction of seasonal abundances.

There is something strange going on with GB in the 2000’s where the LFI increases for a time, but this is no longer the case.

Upper Percentile Drift of the Size distribution

These metrics are an attempt to capture where movements at the edges of the size distribution: Is the size distribution becoming wider or more narrow?

There is a sentiment that the 90th percentile fish is smaller now than it has been in the past. This is true in some circumstances but not all. Indeed, the 90th percentile mass for fishes in the GOM was nearly 2x what it is today. A similar change in the upper percentile location is visible in Georges Bank as well. At the shelf-wide scale and in the southern regions this is less true.

10th percentile weights seem stable and are likely less informative due to gear selectivity preventing consistent sampling at the smallest sizes. I find it more probable that the size distribution changes are related to fluctuations in larger individuals than those at the lower end of the distribution.

# Do them separately
yearly_size_quants <- bind_rows(
    mutate(wigley_in, 
           area = "Northeast Shelf"),
    wigley_in) %>% 
  unite(col = "group_var", area, est_year, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "group_var") %>% 
  separate(group_var, into = c("area", "est_year"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    area = factor(area, levels = area_levels_long_all)) %>% 
  arrange(area)


# Reshape to make filtering easier
quants_long <- yearly_size_quants %>% 
  pivot_longer(
    cols = -c(area, est_year), 
    names_to = "Percentile", 
    values_to = "biomass_kg") 


#-------------

# Take the data excluding dogfish
no_spiny <- filter(wigley_in, comname != "spiny dogfish")


# Do them separately
spineless_size_quants <- bind_rows(
    mutate(no_spiny, area = "Northeast Shelf"),
    no_spiny) %>% 
  unite(col = "group_var", area, est_year, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1,  0.5, 0.90)) %>% 
      t() %>% 
      as_tibble() 
    
  }, .id = "group_var") %>% 
  separate(group_var, into = c("area", "est_year"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    area = factor(area, levels = area_levels_long_all)) %>% 
  arrange(area)


# Reshape to make filtering easier
spineless_quants_long <- spineless_size_quants %>% 
  pivot_longer(
    cols = -c(area, est_year), 
    names_to = "Percentile", 
    values_to = "biomass_kg")

# Separate axes
quants_10 <- quants_long %>% filter(Percentile == "10%") %>% mutate(Percentile = "10th Percentile")
spiny_10 <-  spineless_quants_long %>% filter(Percentile == "10%") %>% mutate(Percentile = "10th Percentile")
q10_p <- ggplot() +
  geom_point(
    data = quants_10,
     aes(est_year, biomass_kg*1000, 
         color = "Wigley Community"),
     size = 0.8, alpha = 0.5) +
  geom_ma(
    data = quants_10,
    aes(x = est_year, y = biomass_kg*1000, 
        color = "Wigley Community",
        linetype = "5-Year Moving Average"),
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  geom_point(
    data = spiny_10,
    aes(est_year, biomass_kg*1000, 
        color = "Spiny Dogfish Removed"),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = spiny_10,
    aes(x = est_year, y = biomass_kg*1000,
        color = "Spiny Dogfish Removed",
        linetype = "5-Year Moving Average"),
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  scale_color_grey() +
  scale_y_continuous(labels = label_number(suffix = "g")) +
  guides(
    linetype = guide_legend(title.position = "top", title.hjust = 0.5),
    color = guide_legend(
      title.position = "top",
      title.hjust = 0.5,
      nrow = 2)) +
  theme(
    strip.text.y = element_blank(),
    plot.margin = margin(0,0,0,0)) +
  facet_grid(area~Percentile, labeller = label_wrap_gen(width = 8)) +
  labs(
    y = "Bodymass", 
    x = "Year", 
    linetype = "Linetype",
    color = "Community")



# Separate axes for 90
quants_90 <- quants_long %>% filter(Percentile == "90%") %>% mutate(Percentile = "90th Percentile")
spiny_90 <-  spineless_quants_long %>% filter(Percentile == "90%") %>% mutate(Percentile = "90th Percentile")
q90_p <-  ggplot() +
  geom_point(
    data = quants_90,
     aes(est_year, biomass_kg, 
         color = "Wigley Community"),
     size = 0.8, alpha = 0.5) +
  geom_ma(
    data = quants_90,
    aes(x = est_year, y = biomass_kg, 
        color = "Wigley Community",
        linetype = "5-Year Moving Average"),
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  geom_point(
    data = spiny_90,
    aes(est_year, biomass_kg, 
        color = "Spiny Dogfish Removed"),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = spiny_90,
    aes(x = est_year, y = biomass_kg,
        color = "Spiny Dogfish Removed",
        linetype = "5-Year Moving Average"),
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  scale_color_grey() +
  scale_y_continuous(labels = label_number(suffix = "kg")) +
  guides(
    linetype = guide_legend(title.position = "top", title.hjust = 0.5),
    color = guide_legend(
      title.position = "top",
      title.hjust = 0.5,
      nrow = 2)) +
  theme(plot.margin = margin(0,0,0,0)) +
  facet_grid(area~Percentile, labeller = label_wrap_gen(width = 8)) +
  labs(
    y = "", 
    x = "Year", 
    linetype = "Linetype",
    color = "Community")



# Assemble
(q10_p | q90_p) + plot_layout(guides = "collect", axes = "collect_x") 

Dogfish & Seasonal Size Distributions

From exploratory analyses we know that spiny dogfish are a large fraction of the biomass in some seasonal communities. Because of this, their behavior has a large impact on the community size structure. However despite their outsized impact they are not the whole story and there is a danger in over-emphasizing their role as they are not strictly benthic and are known to migrate within and beyond the bounds of the survey sampling area.

By excluding spiny dogfish from the community the 2000’s era increase in 90th percentile size for the Gulf of Maine goes away. Possibly giving support to sentiments that spiny dogfish are occupying the body-size-niche perviously held by larger gadids.

# remake the above plot as seasonal

# Do them separately
wigley_season_qs <- bind_rows(
    mutate(wigley_in, 
           area = "Northeast Shelf"),
    wigley_in) %>% 
  unite(col = "group_var", area, est_year, season, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1, 0.5, 0.90)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "group_var") %>% 
  separate(group_var, into = c("area", "est_year", "season"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) %>% 
  arrange(area) %>% 
  pivot_longer(
    cols = -c(area, est_year, season), 
    names_to = "Percentile", 
    values_to = "biomass_kg")  %>% 
  group_by(area, season, Percentile) %>% 
  arrange(est_year) %>% 
  mutate(kg_roll = zoo::rollapply(biomass_kg, 5, mean, na.rm = T, align = "right",  fill = NA))


#-------------

# Take the data excluding dogfish
no_spiny <- filter(wigley_in, comname != "spiny dogfish")


# Do them without dogfish separately
spiny_season_qs <- bind_rows(
    mutate(no_spiny, area = "Northeast Shelf"),
    no_spiny) %>% 
  unite(col = "group_var", area, est_year, season, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1, 0.5, 0.90)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "group_var") %>% 
  separate(group_var, into = c("area", "est_year", "season"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) %>% 
  arrange(area) %>% 
  pivot_longer(
    cols = -c(area, est_year, season), 
    names_to = "Percentile", 
    values_to = "biomass_kg") %>% 
  group_by(area, season, Percentile) %>% 
  arrange(est_year) %>% 
  mutate(kg_roll = zoo::rollapply(biomass_kg, 5, mean, na.rm = T, align = "right",  fill = NA))

Seasonal Size Percentiles

The seasonality in the 90th percentile shift with and w/o spiny dogfish does not speak to migration as much as I thought it might originally.

It is capturing a more general absence of dogfish (and other larger individuals) in MAB during the fall. The departure of larger individuals from the area seasonally is possibly related to the breakdown of the Mid-Atlantic Cold pool, which breaks down over the summer.

# Separate axes for 90
quants_90 <- wigley_season_qs %>% 
  filter(Percentile == "90%") %>% 
  mutate(Percentile = "90th Percentile")
spiny_90 <-  spiny_season_qs %>% 
  filter(Percentile == "90%") %>% 
  mutate(Percentile = "90th Percentile")


# 90th percentile dogfish changes
q90_df_p <-  ggplot() +
  geom_area(
    data = quants_90,
    aes(x = est_year, y = kg_roll, 
        fill = "Complete Community"),
    alpha = 0.6, color = "gray70") +
  geom_area(
    data = spiny_90,
    aes(x = est_year, y = kg_roll, 
        fill = "Spiny Dogfish Removed"),
    alpha = 0.6, color = "gray20") +
  geom_image(
    data = filter(
      quants_90, 
      est_year == 2008, 
      !(area == "Mid-Atlantic Bight" & season == "Fall")), 
    aes(est_year, kg_roll*1.65, 
        image = here::here("Data/phylopics/squalus_phylopic.png")),
    size = 0.35) +
  scale_fill_manual(values = c("gray70", "gray20")) +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  scale_y_continuous(labels = label_number(suffix = "kg")) +
  #coord_flip() +
  guides(
    linetype = guide_legend(title.position = "top", title.hjust = 0.5),
    color = guide_legend(
      title.position = "top",
      title.hjust = 0.5,
      nrow = 2)) +
  theme(plot.margin = margin(0,0,0,0),
        strip.background = element_rect(color = "white", linewidth = 0.2)) +
  facet_nested(
    area~Percentile*season,
    labeller = labeller(
      area = label_wrap_gen(width = 8),
      Percentile = ~str_c(.x, " of Body Mass Distribution"))) +
  labs(
    y = "Individual Bodymass", 
    x = "Year", 
    fill = "Community")
q90_df_p

The above figure shows that the community distribution without dogfish widens in the spring and narrows in the fall.

Percent Dogfish

It is maybe worth noting just how much of the community biomass is held in spiny dogfish, which is frequently above 50%.

dogfish_fraction <- bind_rows(
    mutate(
      wigley_in, 
      area = "Northeast Shelf"),
    wigley_in) %>% 
  group_by(area, est_year, season) %>% 
  summarise(
    dogfish_abund = sum(numlen_adj * (comname == "spiny dogfish"), na.rm = T),
    total_abund   = sum(numlen_adj, na.rm = T),
    dogfish_bio   = sum(sum_weight_kg * (comname == "spiny dogfish"), na.rm = T),
    total_bio     = sum(sum_weight_kg, na.rm = T),
    .groups = "drop") %>% 
  mutate(
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")),
    dogfish_abund_fraction = (dogfish_abund / total_abund),
    dogfish_bio_fraction = (dogfish_bio / total_bio))

dogfish_fraction %>% 
  pivot_longer(cols = ends_with("fraction"), names_to = "metric", values_to = "fraction") %>% 
  mutate(metric = if_else(str_detect(metric, "abund"), "Abundance", "Biomass")) %>% 
ggplot(aes(est_year, fraction)) +
  geom_area(aes(color = season, fill = season), alpha = 0.2) +
  #geom_line(aes(color = season, fill = season), alpha = 0.2) +
  facet_nested(area~metric*season) +
  scale_y_continuous(
    labels = label_percent(),
    limits = c(0,1), expand = expansion(add = c(0,0))) +
  scale_color_gmri() +
  scale_fill_gmri() +
  labs(y = "Community Fraction",
       x = "Year",
       title= " Dogfish as a Fraction of Total Abundance/Bio",
       color = "Season",
       fill = "Season")

Bottom Temperature CCF

Shelf-Wide:

Spring bottom temperatures lead size distributions by 4 or 3 years. The spring temperature part may be important because that is the annual low i.e. the baseline is shifting.

Regionally:

Bottom temperature appears fairly centered (sequential significant correlations at or around 0 lag/lead) in GOM & GB where it is more highly correlated and there is more consistent correlation direction and magnitude at sequential lags.

For MAB-Spring a 1 year of bottom temperature may be important.

In the Fall the GOM has its largest negative correlation centered at a 2-year lead.

# Function to grab the correlation data and lag data
get_ccf_vector <- function(x,y, lagmax = 10){
  
  # Run the ccf
  ccf_data <- ccf(x, y, plot= F , lag.max = lagmax)
  
  # Get the signif:
  # https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/
  # Not 100% sure n is the same for ccf as it is for acf, but...
  significance_level <- qnorm((1 + 0.95)/2)/sqrt(sum(!is.na(x)))
  
  data.frame(
    "acf"    = ccf_data$acf,
    "lag"    = ccf_data$lag,
    "sigpos" = significance_level,
    "signeg" = significance_level*-1
  )
}


# Run Local CCF on Bottom Temperature
btemp_ccf <- wigley_bmspectra_model_df %>% 
  unite(col = "regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(~get_ccf_vector(
    x = .x$bot_temp, 
    y = .x$b), .id = "regseas") %>% 
  separate(col = "regseas", into = c("area", "season"), sep = "XX") %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")))



ggplot(btemp_ccf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(btemp_ccf, sig_flag), color = "black") +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Bottom Temperature CCF",
    x = "Bottom Temperature Lag (years)", 
    y = "ACF", )

Landings CCF

Its plausible that landings have a lagged effect so here are the CCFs for that.

These visually show a pattern of sustained autocorrelation, but in most instances the correlation to landings is such that the slope of the size distribution changes leads the landings AND that higher landings are positively correlated with shallow spectra.

# rejoin other shelf-wide things to shelf-wide landings
shelf_model_df <- wigley_bmspectra_model_df %>% 
  filter(area == "Northeast Shelf") %>% 
  select(-total_weight_lb) %>% 
  left_join(shelf_landings_annual)

# Run Local CCF on Bottom Temperature
landings_ccf <- wigley_bmspectra_model_df %>% 
  filter(area != "Northeast Shelf") %>% 
  bind_rows(shelf_model_df) %>% 
  drop_na(total_weight_lb) %>% 
  unite(col = "regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(~get_ccf_vector(
    x = .x$total_weight_lb, 
    y = .x$b), .id = "regseas") %>% 
  separate(col = "regseas", into = c("area", "season"), sep = "XX") %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")))



ggplot(landings_ccf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(landings_ccf, sig_flag), color = "black") +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Total Landings CCF",
    x = "Total Landings Lag (years)", 
    y = "ACF", )

Spectra Slope ACF

Just as an important check, what is the temporal autocorrelation structure in spectra slopes?

There is a significant positive correlation at 1 year lead in 4/10 regional timeseries. I think the nature of the data is enough to justify adding the AR1() error structure, and this is confirms its inclusion over a year random effect.

There is also six-year lead in 3 of them which is curious, and could be explored as to what is going on there with a 6-year periodicity.

Then lastly the Fall Gulf of Maine spectra has significant correlations up to three years out.

I am of two minds about this:

1. This could be a clue to some state change, as increased autocorrelation is sometimes used in the literature as a signal of hysteresis.

2. OR This could be evidence that Gulf of Maine’s fall community is the exception, and is the ONLY community sampled sufficiently/consistently.

We are in a sense (hopefully) conducting a repeated measure of some community in time, and so there should hopefully be some sort of autocorrelation. The absence of autocorrelation in these communities could

# Pull out acf stuff
spectra_acf <- wigley_bmspectra_model_df %>% 
  unite("regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(function(x){
    # drop NA
    x <- arrange(x, est_year) %>% drop_na(b)
    
    # Do ACF
    x_acf <- acf(x$b, plot = F, lag.max = 10)
    
    # Get the signif:
    # https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/
    # Not 100% sure n is the same for ccf as it is for acf, but...
    significance_level <- qnorm((1 + 0.95)/2)/sqrt(sum(!is.na(x$b)))
    
    data.frame(
      "acf"    = x_acf$acf,
      "lag"    = x_acf$lag,
      "sigpos" = significance_level,
      "signeg" = significance_level*-1
    )
    
  }, .id = "regseas") %>% 
  separate("regseas", into = c("area", "season"), sep = "XX")  %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) 



# Plot things
ggplot(spectra_acf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(spectra_acf, sig_flag), color = "black") +
  scale_x_continuous(
    limits = c(-1,10),
    breaks = c(1:9),
    expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Spectra Slope ACF",
    x = "Lag (years)", 
    y = "ACF", )

Spring and Fall Correlated ??

I can’t believe I never checked this, but are spring and fall communities correlated within an area? My sense is that they should be, and it might be meaningful to point out if they are not.

Only one region (GOM) has “memory” of the community six months before, and has a significant (but small) correlation between spring and fall communities.

Taken with the ACF info, there is consistency that the Gulf of Maine spectra exhibits more sign of temporal autocorrelation and the spectra between season are correlated.

This does however add doubt to how useful these other indices are, in the sense that they may not be capturing change of a community over time.

wigley_bmspectra_model_df %>% 
  select(est_year, area, season, b) %>% 
  pivot_wider(names_from = "season", values_from = "b") %>% 
  ggplot(aes(Spring, Fall)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("P", "R2")), formula = y ~ x) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.5))) +
  facet_wrap(~area, ncol = 2, scales = "free") +
  labs(title = "Spring-Fall Spectra Correlation")

Summary of Trends

Size spectra steepening and median size declines are fairly consistent with one another. Long-term trends are only present at regional scales, and not ocurring shelf-wide for this community of 74 species.

There is something unique happening during the spring in the Mid-Atlantic Bight, where median size is growing and the size spectra is shallowing. This region exhibits large seasonal differences in community composition and large fluctuations in bottom temperature partly due to the mid-atlantic cold-pool.

From the large fish indices I see evidence that individuals within each species are smaller now than they had been, and that larger individuals are rarer in the overall community.

The shift in the 90th percentile size also supports this idea of smaller individuals being more common. The differences seen when excluding spiny dogfish gives some additional evidence that this is particularly true for non-dogfish species. This may also support the idea that dogfish species replacement has ocurred and spiny dogfish are now occupying that larger size class niche space.

From the temperature CCf plot my takeaways are that spring bottom temperature matters more for the shelf-wide size distribution and may operate with a lag of several years. Regionally, the lag seems less relevant and a negative correlation is more evident in GOM+GB.

From the landings CCF I get the sense that total landings aren’t a major driver of slope steepening in our area, or they operate at very long lead times 8-10 years. Many of the significant correlations are in the opposite direction of how we think landings would impact the size distribution (removal of larger individuals leads to steepening), AND they are such that they follow the size distribution slope rather than lead it.

---

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

Temperature and Size Structure

There are known relationships between inter-specific and intra-specific body size across temperature ranges and across latitudinal ranges.

• Bergmann’s Rule: Warmer climates are often inhabited by smaller-bodied species
• James Rule: At an intraspecific level, populations living in warmer conditions often reach smaller maximum body sizes

The metabolic theory of ecology provides a basis for “why” the above rules might be observed.

• Metabolic Theory of Ecology: Metabolic theory predicts how metabolic rate, by setting the rates of resource uptake from the environment and resource allocation to survival, growth, and reproduction, controls ecological processes at all levels of organization from individuals to the biosphere.

Seasonal Movements and Metabolic Efficiency: MTE details how metabolic rates, which are influenced by body size and temperature, are central to the ecological behavior of organisms.

This theory would suggest that larger individuals, which naturally have lower per-unit mass metabolic rates, may move to cooler areas seasonally to optimize their metabolic function and to reduce the energetic costs of thermoregulation in warmer periods.

From what I see at first-glance the size distribution along the shelf is responding consistently with temperature changes (Spring to Fall) on shorter time scales than hypotheses that predict temperature related changes in growth.

wigley_bmspectra_model_df %>% 
  filter(survey_area != "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b, color = season)) +
  geom_point(aes(shape = survey_area)) +
  geom_smooth(method = "lm", color = "black") +
  ggforce::geom_mark_ellipse(aes(fill = season), alpha = 0.1) +
  scale_color_gmri() +
  scale_fill_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  guides(
    shape = guide_legend(nrow = 2),
    fill = guide_legend(nrow = 2),
    color = guide_legend(nrow = 2)) +
  theme(legend.box = "horizontal") +
  labs(
    title = "Bergmann's Rule + MTE",
    subtitle = "The broad relationship between temperature and size distribution is bolstered by\nseasonal community changes in alignment with seasonal temperature differences.",
    shape = "Survey Area",
    fill = "Season",
    color = "Season",
    x = "Bottom Temperature")

Importance of Seasonality

The extent that the regional spectra trends follow a negative linear relationship to increasing temperature is largely due to seasonal differences.

We know that individuals may migrate and/or occupy different habitats seasonally, leading to situations where the community compositions from spring to fall are an apples to oranges comparison. Seasonal movements may be be partially explained by expectations from the MTE, but they also could be separate ecological processes like seasonal foraging or spawning opportunities.

Either way, seasonal differences are large and they justify treating survey seasons as separate groups with independent trends.

Anomalies Model

I wanted to try and avoid conflating that a relationship between spectra and bottom temperature means that “warming” has altered the size spectra. The observation that body size distribution broadly follows some temperature relationship is a confirmation of these broadly observed rules, but I would argue it does not convey that warming has shifted a local climate in some direction.

This is the model I would suggest gets the most at “Warming” impacting regional size distributions:

Fixed effects:

• Season
• Survey area
• Bottom temperature anomaly

Error structure: AR1 autocorrelative errors

By using the season * region interaction groups, we can characterize the typical bottom temperature for each combinations using the long-term seasonal average, and quantify warming using departures from those long-term averages.

The model looks like this as an equation (90% sure), where \(\alpha\) is the intercept, \(\beta x_{tjs}\) is our coefficient for a continuous independent variable (ex. bottom temperature), and \(e_{tjs\)` is our error:

\[y_{tjs} = \alpha + \beta x_{tjs} + e_{tjs}\] and our error depends on the previous timestep:

\[e_{tjs} = \phi e_{t-1,js} + w_{tjs}\] With \(w_{tjs}\) adding noise.

\[w_{tjs} \sim N(0, \sigma^2) \] All of this happening for each \(year_t\) within \(region_j\) & \(season_s\).

I need to consult the Zuur paper about equations in ecology again, but this seems close.

# Just drop the random effect part for now
anom_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp_anom), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(anom_model_ar)$corStruct


# # check the model
# check_model(anom_model_ar)
# plot(check_collinearity(anom_model_ar))
# plot(
#   check_collinearity(
#     nlme::gls(
#       b ~ area + season + scale(bot_temp_anom), #+ scale(log(total_weight_lb)), 
#       data = wigley_bmspectra_model_df, 
#       correlation = corAR1(form = ~ est_year | area/season)
# ))
#   )


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    anom_model_ar, 
    terms = ~ bot_temp_anom * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = anom_model_ar,
  specs =  ~ area * season,
  var = "bot_temp_anom") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp_anom)-2,
    max_temp = max(bot_temp_anom)+2,
    .groups = "drop")



# Plot over observed data
local_btemp_anoms <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp_anom, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp_anoms

Based on this model I would say there is evidence of a warming effect in:
1. GOM-Spring
2. GB-Spring
3. GB-Fall
4. SNE-Fall

Bottom Temperature Model

We were originally using absolute temperatures, which are scaled, and I was curious with the AR error structure how different this would be than the anomaly model.

I just had to see if the same model is better performing or leads to different results when it knows absolute temperatures.

# Just drop the random effect part for now
temp_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(temp_model_ar)$corStruct
# 
# # check the model
# check_model(temp_model_ar)

# # Collinearity without interactions
# plot(
#   check_collinearity(
#      nlme::gls(
#         b ~ area + season + scale(bot_temp), 
#         data = wigley_bmspectra_model_df, 
#         correlation = corAR1(form = ~ est_year | area/season))
# ))


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    temp_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = temp_model_ar,
  specs =  ~ area * season,
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")



# Plot over observed data
local_btemp_results <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Seasonal Bottom Temperature")
local_btemp_results

Verdict: Temperature itself performs better, and is simpler to produce/explain.

data.frame(
  "Model" = c("Temperature Model", "Anomalies Model"),
  "AIC" = c(AIC(temp_model_ar), AIC(anom_model_ar))) %>% 
  mutate(`Delta-AIC` = AIC - min(AIC)) %>% 
  gt()

Model	AIC	Delta-AIC
Temperature Model	-1206.481	0.00000
Anomalies Model	-1184.586	21.89431

Thoughts:

Consistent seasonal shifts in the Shelf-wide spectra suggests to me that the community is more responsive to within-year time scale factors and maybe less sensistive to long-term trends in temperature than I originally expected. This also raises questions about how the spring/fall size distribution differences are being achieved and whether to view them as expressions of the same “community” or separately.

From the decadal variability work we have reason to believe that individual species are both not (perfectly) tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening. These are signals at the species level indicating that temperature effects on growth are ocurring.

In light of that knowledge, there must be some compensatory mechanisms happening across this broader community for the temperature+spectra relationship to be so limited/situational.

Region * Season Pairwise Posthoc

The following figures and tables compare marginal means and post-hoc comparisons for the model above using bottom temperature.

What are the differences across spatial scales and seasonally?

# Region and seasonal posthocs
regseas_phoc <- emmeans(
  temp_model_ar, 
  list(pairwise ~ area + season), 
  adjust = "tukey",
  type = "response")



# Custom Plot
(intercepts_emmeans <- regseas_phoc$`emmeans of area, season` %>% 
    as_tibble() %>%
    ggplot() +
    geom_pointrange(
      aes(area, emmean, ymin = lower.CL, ymax = upper.CL, color = season),
      position = position_dodge(width = 0.25), alpha = 1) +
    #coord_flip() +
    scale_color_gmri() +
    labs(
      y = "Marginal Mean: Mass Spectra Slope",
      x = "Area",
      title = "Wigley Species",
      color = "Season",
      subtitle = "Group marginal means plot"))

# Table for pairwise test results
regseas_phoc %>% 
  pluck(2) %>% 
  as.data.frame() %>% 
  rename(Contrast = `1`) %>% 
  mutate_if(is.numeric, ~round(.x, 3)) %>% 
  arrange(p.value) %>% 
  gt() %>% 
  tab_header(title = "Factor Groups Pairwise Post-Hoc") %>% 
  tab_style(
    style = list(
      cell_fill(color = "#9AFF9A"), 
      cell_text(color = "black")),
    locations = cells_body(
      rows = p.value <= 0.05) 
  )

Factor Groups Pairwise Post-Hoc
Contrast	estimate	SE	df	t.ratio	p.value
Southern New England Spring - Northeast Shelf Fall	0.139	0.030	469.383	4.715	0.000
Southern New England Spring - (Mid-Atlantic Bight Fall)	0.237	0.045	476.537	5.263	0.000
(Mid-Atlantic Bight Spring) - Northeast Shelf Fall	0.157	0.023	455.375	6.921	0.000
(Mid-Atlantic Bight Spring) - (Mid-Atlantic Bight Fall)	0.255	0.041	476.029	6.240	0.000
Northeast Shelf Fall - Gulf of Maine Fall	-0.143	0.029	465.859	-5.004	0.000
Northeast Shelf Fall - Georges Bank Fall	-0.182	0.031	471.728	-5.895	0.000
Gulf of Maine Fall - (Mid-Atlantic Bight Fall)	0.240	0.044	444.388	5.418	0.000
Georges Bank Fall - (Mid-Atlantic Bight Fall)	0.279	0.046	477.818	6.089	0.000
Northeast Shelf Spring - (Mid-Atlantic Bight Fall)	0.201	0.049	473.748	4.104	0.002
Georges Bank Fall - Southern New England Fall	0.121	0.036	477.467	3.412	0.024
Gulf of Maine Spring - (Mid-Atlantic Bight Fall)	0.170	0.050	476.620	3.391	0.026
(Mid-Atlantic Bight Spring) - Southern New England Fall	0.097	0.029	475.519	3.367	0.028
Southern New England Fall - (Mid-Atlantic Bight Fall)	0.158	0.047	415.081	3.335	0.031
Northeast Shelf Spring - Northeast Shelf Fall	0.103	0.035	367.778	2.934	0.100
Georges Bank Spring - (Mid-Atlantic Bight Fall)	0.158	0.057	472.321	2.800	0.139
Gulf of Maine Spring - Georges Bank Fall	-0.109	0.039	477.455	-2.757	0.154
Georges Bank Spring - Georges Bank Fall	-0.121	0.047	476.957	-2.559	0.240
Gulf of Maine Spring - (Mid-Atlantic Bight Spring)	-0.084	0.033	472.244	-2.523	0.259
Gulf of Maine Fall - Southern New England Fall	0.082	0.034	472.581	2.454	0.297
Southern New England Spring - Southern New England Fall	0.079	0.034	476.466	2.294	0.395
Georges Bank Spring - (Mid-Atlantic Bight Spring)	-0.096	0.042	477.515	-2.280	0.404
Northeast Shelf Fall - (Mid-Atlantic Bight Fall)	0.097	0.044	474.070	2.219	0.445
Northeast Shelf Spring - Georges Bank Fall	-0.078	0.038	471.924	-2.073	0.548
Gulf of Maine Spring - Northeast Shelf Fall	0.073	0.037	473.955	1.970	0.621
Gulf of Maine Spring - Gulf of Maine Fall	-0.070	0.038	476.888	-1.853	0.701
Northeast Shelf Fall - Southern New England Fall	-0.060	0.033	471.108	-1.832	0.715
Georges Bank Spring - Gulf of Maine Fall	-0.082	0.046	477.356	-1.790	0.741
Gulf of Maine Spring - Southern New England Spring	-0.066	0.038	477.822	-1.729	0.779
Northeast Shelf Spring - (Mid-Atlantic Bight Spring)	-0.054	0.031	464.727	-1.716	0.786
Georges Bank Spring - Southern New England Spring	-0.078	0.046	477.143	-1.694	0.799
Georges Bank Spring - Northeast Shelf Fall	0.061	0.045	474.606	1.348	0.942
Southern New England Spring - Georges Bank Fall	-0.042	0.032	476.714	-1.301	0.953
Gulf of Maine Fall - Georges Bank Fall	-0.039	0.031	440.843	-1.235	0.966
Northeast Shelf Spring - Gulf of Maine Fall	-0.039	0.036	468.487	-1.095	0.985
Northeast Shelf Spring - Southern New England Fall	0.043	0.039	471.482	1.092	0.985
Northeast Shelf Spring - Southern New England Spring	-0.036	0.037	470.472	-0.980	0.993
(Mid-Atlantic Bight Spring) - Georges Bank Fall	-0.024	0.026	477.532	-0.918	0.996
Northeast Shelf Spring - Georges Bank Spring	0.043	0.050	474.220	0.848	0.998
Northeast Shelf Spring - Gulf of Maine Spring	0.030	0.043	473.560	0.709	0.999
Southern New England Spring - (Mid-Atlantic Bight Spring)	-0.018	0.025	434.985	-0.723	0.999
Gulf of Maine Spring - Georges Bank Spring	0.012	0.051	422.378	0.235	1.000
Gulf of Maine Spring - Southern New England Fall	0.013	0.041	474.664	0.308	1.000
Georges Bank Spring - Southern New England Fall	0.001	0.049	437.668	0.012	1.000
Southern New England Spring - Gulf of Maine Fall	-0.003	0.030	477.749	-0.109	1.000
(Mid-Atlantic Bight Spring) - Gulf of Maine Fall	0.015	0.024	472.981	0.622	1.000

regseas_phoc %>% 
  pluck(2) %>% 
  as.data.frame() %>% 
  rename(Contrast = `1`) %>% 
  mutate_if(is.numeric, ~round(.x, 3)) %>% 
  filter(p.value < 0.05)  %>% 
  separate(Contrast, into = c("group1", "group2"), sep = " - ") %>% 
  mutate(
    group1 = str_remove_all(group1, "[(]|[)]"),
    group2 = str_remove_all(group2, "[(]|[)]"),
    statement = 
      if_else(estimate > 0, " is shallower than ", " is steeper than ")) %>% 
  select(group1, statement, group2) %>% 
  arrange(statement, group1) %>% 
  gt()

group1	statement	group2
Georges Bank Fall	is shallower than	Southern New England Fall
Georges Bank Fall	is shallower than	Mid-Atlantic Bight Fall
Gulf of Maine Fall	is shallower than	Mid-Atlantic Bight Fall
Gulf of Maine Spring	is shallower than	Mid-Atlantic Bight Fall
Mid-Atlantic Bight Spring	is shallower than	Northeast Shelf Fall
Mid-Atlantic Bight Spring	is shallower than	Southern New England Fall
Mid-Atlantic Bight Spring	is shallower than	Mid-Atlantic Bight Fall
Northeast Shelf Spring	is shallower than	Mid-Atlantic Bight Fall
Southern New England Fall	is shallower than	Mid-Atlantic Bight Fall
Southern New England Spring	is shallower than	Northeast Shelf Fall
Southern New England Spring	is shallower than	Mid-Atlantic Bight Fall
Northeast Shelf Fall	is steeper than	Gulf of Maine Fall
Northeast Shelf Fall	is steeper than	Georges Bank Fall

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

# Plot annual totals
bind_rows(list(
  shelf_landings_annual,
  landings_annual)) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_long_all)) %>% 
  filter(est_year <= 2019) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.8) +
  geom_ma(
    aes(color = "5-Year Moving Average"), 
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries", linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1, scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Landings Autoregressive model

I wanted to follow on the bottom temperature anomalies model here for landings:

# Just drop the random effect part for now
f_model_ar <- nlme::gls(
  b ~ area * season + area * scale(log(total_weight_lb)), 
  data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>% 
    dplyr::filter(est_year > 1978),
  correlation = corAR1(form = ~ est_year | area/season)
)



# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    f_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = f_model_ar,
  ~area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite")


# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = f_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(mod2_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free") +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

# temp model is better
# AIC(temp_model_ar)
# AIC(f_model_ar)

This is our old friend “higher landings when spectra was shallower” that we’ve seen before.

Temperature and Landings Together

If we want to just throw them both in because we’ve done what we can and want to hold them both in the frame at once we get this:

# Just drop the random effect part for now
full_model_ar <- nlme::gls(
  # Model
  model = b ~ area * season * scale(bot_temp) + area * scale(log(total_weight_lb)), 
  
  # Data
  data = wigley_bmspectra_model_df %>%
    drop_na(total_weight_lb) %>% 
    dplyr::filter(area != "Northeast Shelf"),
  
  # Autocorrelation
  correlation = corAR1(form = ~ est_year | area/season))

# # Collinearity without interactions
# plot(
#   check_collinearity(
#     nlme::gls(
#         b ~ area + season + scale(bot_temp) +  scale(log_land),
#         data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>%
#                  filter(area != "Northeast Shelf") %>% 
#           mutate(log_land = log(total_weight_lb)),
#         correlation = corAR1(form = ~ est_year | survey_area/season)
#       )
#     )
#   )


# # confidence intervals on phi
# intervals(full_model_ar)$corStruct


# temp model is much better
# AIC(temp_model_ar)
# AIC(full_model_ar)

Before proceeding it is worth noting that this model is less parsimonious than a temperature or anomaly only model:

data.frame(
  "Model" = c(
    "Temperature Model", 
    "Anomalies Model", 
    "Landings Model", 
    "Landings and Temperature"), 
  "AIC" = c(
    AIC(temp_model_ar), 
    AIC(anom_model_ar),
    AIC(f_model_ar),
    AIC(full_model_ar))) %>% 
  mutate(`Delta-AIC` = AIC - min(AIC)) %>% 
  gt()

Model	AIC	Delta-AIC
Temperature Model	-1206.4805	0.00000
Anomalies Model	-1184.5862	21.89431
Landings Model	-744.1868	462.29373
Landings and Temperature	-897.3938	309.08670

# Clean up the plot:
# Plot the predictions over data
temp_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
temp_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area * season,
  mode = "appx-satterthwaite",
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  filter(area != "Northeast Shelf") %>% 
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")



# Plot over observed data
local_btemp <- temp_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(temp_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area!= "Northeast Shelf"),
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp

# Clean up the plot:
# Plot the predictions over data
f_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
f_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite"
  ) %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- f_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(f_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area != "Northeast Shelf"),
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

https://noaa-edab.github.io/tech-doc/zoo_abundance_anom.html

zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo|oith")) %>% 
  mutate(EPU = factor(EPU, levels = c("GOM", "GB", "MAB"))) %>% 
ggplot(aes(Time, Value)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Moving Average"),
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  #scale_color_gmri() +
  facet_grid(EPU~Var) +
  labs()

Cold Pool Index

There are several cold pool indices available that might be helpful in adding context to changes in MAB. Below are the Dupontavice et al. 2022 methods that are used in the SOE report.

The Cold Pool Index (Model_CPI) was adapted from Miller et al. (2016). Residual temperature was calculated in each grid cell, \(i\), in the Cold Pool domain as the difference between the average bottom temperature at the year \(y\) (\(T_{i,y}\)) and the average bottom temperature over the period 1972–2019 (\(\bar{T}_{i, 1972-2019}\)⁠⁠) between June and September. \(CPI\) was calculated as the mean residual temperature over the Cold Pool domain such that

\[CPI_y~=~\frac{\sum_{i=1}^{n}~(T_{i,y}-\bar{T}_{i,~1972-2019})}{n}\]

where n is the number of grid cells over the Cold Pool domain.

Takeaway:

\(CPI\) was calculated as the mean residual temperature over the Cold Pool domain

From this index we can see that the temperatures over the cold pool domain (which include the MAB) have increased over time.

cpool <- ecodata::cold_pool
cpool %>% 
  filter(!str_detect(Var, "se_")) %>% 
  ggplot(aes(Time, Value)) +
  geom_line() +
  facet_grid(Var~EPU, scales = "free")

# # There's an SF
# cpool_sf <- ecodata::cold_pool_sf %>% 
#   mutate(id = row_number())
# 
# ggplot(cpool_sf) + geom_sf(aes(fill = id), alpha = 0.2)

---

Driver Correlations